Localized atomic basis set in the projector augmented wave method 
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We present an implementation of localized atomic orbital basis sets in the projector augmented 
wave (PAW) formalism within the density functional theory (DFT). The implementation in the real- 
space GPAW code provides a complementary basis set to the accurate but computationally more 
demanding grid representation. The possibility to switch seamlessly between the two representations 
implies that simulations employing the local basis can be fine tuned at the end of the calculation 
by switching to the grid, thereby combining the strength of the two representations for optimal 
performance. The implementation is tested by calculating atomization energies and equilibrium 
bulk properties of a variety of molecules and solids, comparing to the grid results. Finally, it is 
demonstrated how a grid-quality structure optimization can be performed with significantly reduced 
computational effort by switching between the grid and basis representations. 

PACS numbers: 71.15.Ap, 71.15.Dx, 71.15.Nc 
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I. INTRODUCTION 

Density functional theory (DFT) with the single- 
particle Kohn-Sham (KS) scheme is presently the most 
widely used method for electronic structure calculations 
in both solid state physics and quantum chemistry^— 
Its success is mainly due to a unique balance between ac- 
curacy and efficiency which makes it possible to handle 
systems containing hundreds of atoms on a single CPU 
with almost chemical accuracy. 

At the fundamental level the only approximation of 
DFT is the exchange-correlation functional which con- 
tains the non-trivial parts of the kinetic and electron- 
electron interaction energies. However, given an 
exchange-correlation functional one is still left with the 
non-trivial numerical task of solving the Kohn-Sham 
equations. The main challenge comes from the very rapid 
oscillations of the valence electrons in the vicinity of the 
atom cores that makes it very costly to represent this 
part of the wavefunctions numerically. In most mod- 
ern DFT codes the problem is circumvented by the use 
of pseudopotentials4~— The pseudopotential approxima- 
tion is in principle uncontrolled and is in general subject 
to transferability errors. An alternative method is the 
projector augmented wave (PAW) method invented by 
BlochI -. An appealing feature of the PAW method is that 
it becomes exact if sufficiently many projector functions 
are used. In another limit the PAW method becomes 
equivalent to the ultra-soft pseudopotentials introduced 
by Vanderbilt 5 -. 

The representation of the Kohn-Sham wavefunctions 
is a central aspect of the numerics of DFT. High ac- 
curacy is achieved by using system independent ba- 
sis sets such as plane waves 7 . - — , wavelet s 10 ! 11 or real- 
space grid a 12 ' 13 , which can be systematically expanded to 
achieve convergence. Less accurate but computationally 
more manageable methods expand the wavefunction in 
terms of a system-dependent localized basis consisting of 
e.g. Gaussians— or numerical atomic orbital o 15 ' 16 . Such 



basis sets cannot be systematically enlarged in a simple 
way, and consequently any calculated quantity will be 
subject to basis set errors. For this reason the former 
methods are often used to obtain binding energies where 
accuracy is crucial, while the latter are useful for struc- 
tural properties which are typically less sensitive to the 
quality of the wavefunctions. 

In this paper we discuss the implementation of a local- 
ized atomic basis set in the PAW formalism and present 
results for molecular atomization energies, bulk proper- 
ties, and structural relaxations. The localized basis set, 
which we shall refer to as the LCAO basis, is similar to 
that of the well-known Siesta pseudopotential codei£, but 
here it is implemented in our recently developed multi- 
grid PAW code GPAW—. A unique feature of the re- 
sulting scheme is the possibility of using two different 
but complementary basis sets: On the one hand wave- 
functions can be represented on a real-space grid which 
in principle facilitates an exact representation, and on 
the other hand the wavefunctions can be represented in 
the efficient LCAO basis. This allows the user to switch 
seamlessly between the two representations at any point 
of a calculation. As a particularly powerful application 
of this "double-basis" feature, we demonstrate how ac- 
curate structural relaxations can be performed by first 
relaxing with the atomic basis set and then switching to 
the grid for the last part. Also adsorption energies, which 
are typically not very good in LCAO, can be obtained on 
the grid at the end of a relaxation. 

While LCAO pseudopotential codes as well as plane- 
wave/grid PAW codes already exist and have been dis- 
cussed extensively in the literature ^ 15 ' 16 the combina- 
tion of LCAO and PAW is new. Compared to the popular 
Siesta method, which is based on norm-conserving pseu- 
dopotentials, the advantage of the present scheme (apart 
from the "double basis" feature) is that PAW works with 
coarser grids to represent the density and effective po- 
tentials. As an example, Fig.[T]shows the atomic orbitals 
of iron calculated with the norm-conserving Hartwigsen- 
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FIG. 1: (Color online) The pseudo valence states of iron calcu- 
lated with PAW and the norm-conserving HGH pseudopoten- 
tials. Both methods produce smooth wave functions for the 
delocalized 4s state, but the lack of norm conservation allows 
the short-ranged 3d state in PAW to be accurately sampled 
on a much coarser grid. 

Goedecker-Hutter (HGH) pseudopotentials^ as well as 
with PAW. Clearly the d wavefunction is much smoother 
in PAW. This is essential for larger systems where op- 
erations on the grid, i.e. solving the Poisson equation, 
evaluating the density, and calculating the potential ma- 
trix elements, become computationally demanding. 

II. PROJECTOR AUGMENTED WAVE 
METHOD 

In this section we give a brief review of the PAW for- 
malism. For simplicity we restrict the equations to the 
case of spin-paired, finite systems, but the generaliza- 
tions to magnetic and periodic systems are straightfor- 
ward. For a more comprehensive presentation we refer 
to Ref. 0. 

A. PAW Transformation Operator 

The PAW method is based on a linear transforma- 
tion T which maps some computationally convenient 
"pseudo" or "smooth" wavefunctions \ip n ) to the phys- 
ically relevant "all-electron" wavefunctions \ip n )- 

\^n)=T\i> n }, (1) 

where n is a quantum state label, consisting of a band 
index and possibly a spin and k-vector index. 

The transformation is chosen as T = T a , i.e. the 

identity operator plus an additive contribution centered 
around each atom, which differs based on the species of 
atom. 

The atomic contribution for atom a is determined by 
choosing a set of smooth functions <f>f(r), called pseudo 
partial waves, and requiring the transformation to map 
those onto the atomic valence orbitals <f>f(r) of that 



atom, called all-electron partial waves. This effectively 
allows the all-electron behaviour to be incorporated by 
the smooth pseudo wave functions. Since the all-electron 
wave functions are smooth sufficiently far from the atoms, 
we may require the pseudo partial waves to match the all- 
electron ones outside a certain cutoff radius, such that 
<j>f(r) = <pf(r) for r > r c . This localizes the atomic con- 
tribution T a to the augmentation sphere r < r c . Finally 
a set of localized projectors pf(r) is chosen as a dual basis 
to the pseudo partial waves. We further want the partial 
wave-projector basis to be complete within the augmen- 
tation sphere, in the sense that any pseudo wave function 
should be expressible in terms of pseudo partial waves, 
and therefore require 

£|#><p?| = l, (ft\%)=6 ij . (2) 

i 

The transformation T is then defined by 

r=i+j2H(\ ( i ) t)-\kmi (3) 

a i 

which allows the all-electron Kohn-Sham wavefunction 
V'n(r) = (r|i/ , n) to be recovered from a pseudo wave func- 
tion through 

^(r) = Mr) + ££>?(r) - #(r))<p?$ B ). (4) 

a i 

We emphasize that the all-electron wave functions are 
never evaluated explicitly, but all-electron values of ob- 
servables are calculated through manipulations which 
rely only on coarse grids or one-dimensional radial grids. 
Using using Eqs. ([I]) and (j3|), the all-electron expectation 
value for any semi-local operator O due to the valence 
states can be written 

(o)=Y,umo\i> n ) 

n 

naij 
naij 

Inside the augmentation spheres the partial wave expan- 
sion is ideally complete, so the first and third terms will 
cancel and leave only the all-electron contribution. Out- 
side the augmentation spheres the pseudo partial waves 
are identical to the all-electron ones, so the two atomic 
terms cancel. The atomic matrix elements of O in the 
second and third terms can be pre-evaluated for the iso- 
lated atom on high-resolution radial grids, so operations 
on smooth quantities, like (ip n \0\ip n ) and (pf\ip n ), are 
the only ones performed during actual calculations. 
It is convenient to define the atomic density matrices 

D^=J2(Pt\^n)fn^nm, (6) 
n 
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since these completely describe the dependence of the 
atomic terms in Eq. ([5]) on the pseudo wave functions. 
The expectation value can then be written 

(0)=^f n $ n \0$ n ) 

71 

+E D 5i(wW)-(*?W>)' ( ? ) 

aij 

Although the PAW method is an exact implementation 
of density functional theory, some approximations are 
needed for realistic calculations. The frozen-core approx- 
imation assumes that the core states are localized within 
the augmentation spheres and that they are not modi- 
fied by the chemical environment and hence taken from 
atomic reference calculations. The non-completeness of 
the basis, or equivalently the finite grid-spacing, will in- 
troduce an error in the evaluation of the PS contribu- 
tion ip n in ([5]). Finally, the number of partial waves and 
projector functions is obviously finite. This means that 
the completeness conditions of Eq. ^ we have required 
are not strictly fulfilled. This approximation can be con- 
trolled directly by increasing the number of partial waves 
and projectors. 

B. Density 

The electron density n(r) is the expectation value of 
the real-space projection operator and, by Eq. ([7]), takes 
the form 

n(r) = n(r) + E [n a (r - R a ) - n a {v - R Q )] , (8) 



where 



Here we have separated out the all-electron core density 
n a c (r) and the pseudo core density ri° (r) , where the latter 
can be chosen as any smooth continuation of n"(r) inside 
the augmentation spheres, since it will cancel out in Eq. 
([8]). We omit conjugation of the partial waves since these 
can be chosen as real functions without loss of generality. 



C. Compensation charges 

In order to avoid dealing with the cumbersome nu- 
clear point charges, and to compensate for the lack of 
norm-conservation, we introduce smooth localized com- 
pensation charges Z a (v) on each atom, which are added 



to n(r) and n a (r), thus keeping the total charge neutral. 
This yields a total charge density that can be expressed 
as 

p{v) = p(r) + ]T [p a (r - R a ) - p a (r - R a )] , (12) 

a 

in terms of the neutral charge densities 

p(r) = n(r) + Z(r) = n(r) + E Z a (r - R a ), (13) 



p a (r) = n a (r) + Z a d(r), 
p a (r)=h a (r) + Z a (r), 



(14) 
(15) 



where Z a S(r) is the central nuclear point charge. The 
compensation charges are chosen to be localized functions 
around each atom of the form 

Z a (r) = £ Ql~gl(r) = E Qf m r l gf (r)Y lm (r) , (16) 



where g^(r) are fixed Gaussians, and Yi m (r) are spheri- 
cal harmonics. We use L = I, m as a composite index for 
angular and magnetic quantum numbers. The expansion 
coefficients Q a h are determined in terms of Dfj by requir- 
ing the compensation charges to cancel all the multipolc 
moments of each augmentation region up to some order, 
generally Z max = 2. The charges will therefore dynami- 
cally adapt to the surroundings of the atom. For more 
details we refer to the original work by Blochl 7 . 



D. Total Energy 

The total energy can also be separated into smooth 
and atom-centered contributions 



n(r) = ^/„|^(r)| 2 +^^(|r-R a |), (9) 

n a 

i a (r) = X>^r)^(r)+<(r), (10) where 

ij 

i a (r) = J2D° i 4>?(r)%(r)+n a c (r). (11) 



E = E + J2( E<1 - E<1 ), 



(17) 



^ = E/«^«|-I v2 |^)+E / "(rK(|r-R a |)dr 

n a J 

1 ^^W + lUnL (18) 



2 77 |r - r' 



(19) 



The terms T c a oro and T£ OTe are the kinetic energy con- 
tributions from the frozen core states, while v a (r) is an 
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arbitrary potential, vanishing for r > r". This potential 
is generally chosen to make the atomic potential smooth, 
while its contribution to the total energy vanishes if the 
partial wave expansion is completed 

E xc is the exchange-correlation functional, which must 
be local or semilocal as per Eq. for the above expres- 
sions to be correct. While the functional is non-linear, it 
remains true that 



E xc [n] = E xc [h] 



{E xc [n a } - E xc [fi a }) , (21) 



III. LOCALIZED BASIS SETS IN PAW 

We now introduce a set of basis functions which 
are fixed, strictly localized atomic orbital-like func- 
tions represented numerically, following the approach by 
Sankey and Niklewskii^. We furthermore consider the 
pseudo wave functions \ipn) to be linear combinations of 
the new basis functions 



(28) 



because of the functional's semilocality: the energy con- 
tribution from n(r) around every point inside the aug- 
mentation sphere is exactly cancelled by that of n a (r), 
since n(r) and n a (r) are exactly identical here, leaving 
only the contribution £" xc [n a ]. Outside the augmenta- 
tion region, a similar argument applies to n a (r) and n(r), 
leaving only the energy contribution from n(r) which is 
here equal to the all-electron density. 



E. Hamiltonian and orthogonality 

In generic operator form, the Hamiltonian correspond- 
ing to the total energy from Eq. (IT71) is 



J2\Pt)^H^\, (22) 



where v = «H a [p] + v + w xc [n] is the local effective poten- 
tial, containing the Hartree, the arbitrary localized and 
the xc potentials, and where 



AH a = ■ 

13 dD% 



(23) 



are the atomic Hamiltonians containing the atom- 
centered contributions from the augmentation spheres. 
Since the all-electron wave functions tpn must be or- 
thonormal, the pseudo wave functions ^ n must obey 

Sum = (lpn\lpm) = $n\PT\i> m ) = $n\S$ m ), (24) 
where we have defined the overlap operator 

S = t t T=l + '£ f fflAS? j (fi\. (25) 



The atomic contributions 



AS?- = (tf |# 



(26) 



are constant for a given element. 

Given the Hamiltonian and orthogonality condition, a 
variational problem can be derived for the pseudo wave 
functions. This problem is equivalent to the generalized 
Kohn-Sham eigenvalue problem 



H\tp n ) = S\tp n )e ni 



(27) 



which can then be solved self-consistently with available 
techniques. 



where the coefficients c^ n are variational parameters. It 
proves useful to define the density matrix 



v — ^ CfmfnC^r 



(29) 



The pseudo density can be evaluated from the density 
matrix through 

f»(r) = 2$*(r)* tf (rK M + 5^n»(r). (30) 

Ahead of a calculation, we evaluate the matrices 

= <<M - iV 2 |$,), (31) 



P?» = (p?I*m). 



(32) 
(33) 



which are used to evaluate most of the quantities of the 
previous sections in matrix form. The atomic density 
matrices from Eq. ^) become 



D a — \ " P a n P a * 



(34) 



and the kinetic energy contribution in the first term of 
Eq. CH]) is 



(35) 



We can then define the Hamiltonian matrix elements by 
taking the derivative of the total energy E with respect 
to the density matrix elements, which eventually results 
in the discretized Hamiltonian 



H 



ft V 



dE 



V + V^ + J2 P^&Hf^, (36) 



where 



V IW = J $;(r)5(r)$ v (r)dr. 



(37) 



The overlap operator of Eq. (|25[) has the matrix repre- 
sentation 

S, v = (%\S\$„) = G M „ + PiZ&StjP^, (38) 
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so orthogonality of the wave functions is now expressed 
by 

This is incorporated by defining a quantity f2 to be varia- 
tionally minimized with respect to the coefficients, specif- 
ically 

= E — y ' A nm (c^SfiuCpn — <^mn) • (40) 
mil/if 

Setting the derivative of 17 with respect to c^n equal to 
0, one obtains the generalized eigenvalue equation 

^ ^ H pi/Ci/fi — ^ ' ^i/^n^ni (41) 

which can be solved for the coefficients c^ n and energies 
e n when the Hamiltonian H^ u and the overlap matrix 
S^ u are known. 

A. Basis functions generation 

The basis functions in Eq. (|28p are atom-centered 
orbitals written as products of numerical radial functions 
and spherical harmonics: 

$«im(r) = ip n i(r)Y lm (r). (42) 

In order to make the Hamiltonian and overlap matrices 
sparse in the basis-set representation, we use strictly lo- 
calized radial functions, i.e. orbitals that are identically 
zero beyond a given radius, as proposed by Sankey and 
Niklewski 15 and successfully implemented in the SIESTA 
method^. 

The first (single-zeta) basis orbitals fnF( r ) are 
tained for each valence state by solving the radial all- 
electron Kohn-Sham equations for the isolated atom in 
the presence of a confining potential with a certain cut- 
off. If the confining potential is chosen to be smooth, 
the basis functions similarly become smooth. We use the 
same confining potential as proposed in Ref. UtI The 
smooth basis functions are then obtained using (p„i (r) = 
T _1 <p^; E (r). The result of the procedure is illustrated in 

Fig. m 

The cutoff radius is selected in a systematic way by 
specifying the energy shift AE of the confined orbital 
compared to the free-atom orbital. In this approach 
small values of AE will correspond to long-ranged ba- 
sis orbitals^. 

To improve the radial flexibility, extra basis functions 
with the same angular momentum I (multiple-zeta) are 
constructed for each valence state using the split- valence 
technique^. The extra function is constructed by match- 
ing a polynomial to the tail of the atomic orbital, where 
the matching radius is determined by requiring the norm 
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FIG. 2: NAO generation for the nitrogen 2s state: the all- 
electron orbital of the free atom, the confined all-electron or- 
bital, and the corresponding pseudo wave function after ap- 
plying the inverse PAW transformation. The augmentation 
sphere and basis function cutoffs are indicated. 



of the part of the atomic orbital outside that radius to 
have a certain value. 

Finally, polarization functions (basis functions with I 
quantum number corresponding to the lowest unoccupied 
angular momentum) can be added in order to improve 
the angular flexibility of the basis. There are several 
approaches to generate these orbitals, such as perturb- 
ing the occupied eigenstate with the highest / quantum 
number with an electric field using first order pertur- 
bation theory (like in Ref. HH) or using the appropriate 
unoccupied orbitals. As a first implementation we use a 
Gaussian-like function of the form r l exp(—ar 2 ) for the 
radial part, where I corresponds to the lowest unoccupied 
angular momentum. This produces reasonable polariza- 
tion functions as demonstrated by the results presented 
in a following section. 

A generator program is included in the GPAW code 
and it can produce basis sets for virtually any elements 
in the periodic table. Through our experiences with gen- 
erating and using different basis sets, we have reached 
the following set of default parameters: We usually work 
with a DZP basis. The energy shift for the atomic orbital 
is taken as 0.1 eV, and the tail-norm is 0.16 (in agree- 
ment with SIESTA^). The width of the Gaussian used 
for the polarization function is 1 /4 of the cut-off radius of 
the first zeta basis function. Further information can be 
found in the documentation for the basis set generator. 
At this point we have not yet systematically optimized 
the basis set parameters, although we expect to do so by 
means of an automatic procedure. 



B. Atomic forces 

The force on some atom a is defined as the negative 
derivative of the total energy of the system with respect 
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to the position of that atom, 



F a = - 



dE 



(43) 



The derivative is to be taken with the constraints that 
selfconsistency and orthonormality according to (|31?)) 
must be obeyed. This implies that the calculated force 
will correspond to the small-displacement limit of the 
finite-difference energy gradient one would obtain by per- 
forming two separate energy calculations, where atom a 
is slightly displaced in one of them. 

The expression for the force is obtained by using the 
chain rule on the total energy of Eq. (1171) . The primary 
complication compared to the grid-based PAW force for- 
mula, Eq. (50) from Ref. [H, is that the basis functions 
move with the atoms, introducing extra terms in the 
derivative. 

The complete formula for the force on atom a is 

F " = m E Sta* - ^ E 



25ft 22, ZltuEvt* - 2jfty^Z^E, 



UfX 



fJLU 



2S £.1/ 



d$* (r) 
^z)(r)^(r)dr 



dn"(|r 



v(r)- 



dR' 

-R a l) 



Pun 



dR a 



dr 



^(|r-R a |) 
n(r) r^rr; dr 



dR a 



L 



dR a 



where 



A 6 - 



r] P b * 



dR a 



E 



dP b * 

__lfJ_ A erf) pb 



— 



(44) 

(45) 

(46) 
(47) 



The notation \i S a denotes that summation should be 
performed only over those basis functions that reside on 
atom a. 

Eq. (|4"4")l is derived in Appendix|XJ The last three terms 
are basis set independent, and inherited from the grid- 
based implementation. 



IV. IMPLEMENTATION 

The LCAO code is implemented in GPAW, a real-space 
PAW code. For the details of the real-space implemen- 
tation we refer to the original paper—. In this code the 



density, effective potential and wave functions are evalu- 
ated on real-space grids. 

In LCAO the matrix elements of the kinetic and over- 
lap operators T^, 6 M „ and P? in Eqs. (|3"Tj) - (|3"3"]) are 
efficiently calculated in Fourier space based on analytical 
expressions^. For each pair of different basis orbitals 
(i.e. independently of the atomic positions), the over- 
lap can be represented in the form of radial functions 
and spherical harmonics. These functions are stored as 
splines which can in turn be evaluated for a multitude of 
different atomic separations. 

The two-center integrals are thus calculated once for a 
given atomic configuration ahead of the self-consistency 
loop. This is equivalent to the SIESTA approach^. 

The matrix elements of the effective potential V^ v are 
still calculated numerically on the three dimensional real- 
space grid, since the density is also evaluated on this 
gricU^. 

Because of the reduced degrees of freedom of a ba- 
sis calculation compared to a grid-based calculation, the 
Hamiltonian from Eq. (|36[) is directly diagonalized in 
the space of the basis functions according to Eq. (|41l) . 
This considerably lowers the number of required itera- 
tions to reach selfconsistency, compared to the iterative 
minimization schemes used in grid-based calculations. 

For each step in the selfconsistency loop, the Hartree 
potential 5n a (r) is calculated by solving the Poisson 
equation V 2 VHa(r) = — 47r j 5(r) in real space using existing 
multigrid methods, such as the Gauss-Seidel and Jacobi 
methods. A solver based on the fast Fourier transform is 
also available in the GPAW code. 

The calculations are parallelized over k-points, spins 
and real-space domains like in the grid-based case^. We 
further distribute the orbital-by-orbital matrices such as 
and S^, and use ScaLAPACK for operations on 
these, notably the diagonalization of Eq. (|4T|) . 



A. Localized functions on the grid 

Quantities such as the density n(r) and effective po- 
tential u(r) are still stored on 3D grids. Matrix elements 
like Vpv in Eq. (I3T1) . and the pseudo density given by 
Eq. (|30)) , can therefore be calculated by loops over grid 
points. 

Since each basis function is nonzero only in a small 
part of space, we only store the values of a given function 
within its bounding sphere. Each function value inside 
the bounding sphere is calculated as the product of ra- 
dial and angular parts vz. Eq. (|42[) . where the radial part 
is represented by a spline, and the spherical harmonic 
evaluated in cartesian form, i.e. as a polynomial. The 
same method is used to evaluate derivatives in force cal- 
culations, although this involves the derivatives of these 
quantities aside from just their function values. 

We initially compile a data structure to keep track of 
which functions are nonzero for each grid point. When 
looping over the grid, we maintain a list of indices \x for 
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the currently nonzero functions by adding or removing, 
as appropriate, those functions whose bounding spheres 
we intersect. The locations of these bounding spheres 
are likewise precompiled into lists for efficient process- 
ing. The memory overhead due to this method is still 
much smaller than the storage requirements for the ac- 
tual function values. 



V. RESULTS 

In this section we calculate common quantities using 
the localized basis set on different systems. The results 
are compared to the complete basis set limit, i.e. a well 
converged grid calculation. Note that this comparison 
can be done in a very systematic way since the calcu- 
lations on the grid share the same approximations and 
mostly the same implementation as the calculations per- 
formed with the localized basis. All the results presented 
in this section have been obtained using PAW setups from 
the extensive GPAW library, freely available online^. 



A. Molecules 

In order to assess the accuracy of the LCAO imple- 
mentation for small molecules, the PBE^i. atomization 
energies for the G2-1 data-set^ are considered. The 
atomic coordinates are taken from MP2(full)/6-31G(d) 
optimized geometries. The error with respect to the grid 
results is shown in Fig. [3] for different basis sets. This 
error is defined as 

AE LCAO AES ri d = £LCAO _ £ ^LCAO 

atoms 

- - E ( 48 ) 

\ atoms / 

The reference grid results are well converged calculations 
in very good agreement with the VASP 8 and Gaussian 14 
codes. 

The figure shows that enlarging the basis set, i.e. in- 
cluding more orbitals per valence electron, systematically 
improves the results towards the grid energies. 

It must be noted that some differences with respect 
to the grid atomization energies still remain, even in the 
case of large basis sets. This is mainly due to the two fol- 
lowing reasons. Firstly, the basis functions are generated 
from spin-paired calculations and hence they do not ex- 
plicitly account for possible spin-polarized orbitals. This 
is in practice accounted for by using larger basis-sets in 
order to include more degrees of freedom in the shape 
of the wavefunctions. Secondly, isolated atoms are diffi- 
cult to treat because of their long ranged orbitals. Actual 
basis functions are, in fact, obtained from atomic calcula- 
tions with an artificial confining potential thus resulting 
in more confined orbitals. 



B. Solids 

The equilibrium bulk properties have been calculated 
for several crystals featuring different electronic struc- 
tures: simple metals (Li, Na, Al), semiconductors (A1P, 
Si, SiC), ionic solids (NaCl, LiF, MgO) transition metals 
(Fe, Cu, Pt) as well as one insulator (C). The results are 
shown in Fig. 2J For comparison with grid-based calcula- 
tions, the bar plots show the deviations from grid-based 
results for each basis set, while the precise numbers are 
shown in each of the corresponding tables. All the cal- 
culations were performed with the solids in their lowest 
energy crystal structure, using the PBE functional for ex- 
change and correlational. The quantities were computed 
using the relaxed structures obtained with the default, 
unoptimized basis sets. The calculations were generally 
spin-paired, i.e. non-magnetic, with the exception of Fe 
and the atomic calculations used to get cohesive energies. 

The overall agreement with the real-space grid is ex- 
cellent: about 0.5% mean absolute error in the compu- 
tation of lattice constants, 4% in cohesive energies and 
5-8% for bulk moduli using DZP basis sets. Notice that in 
many cases remarkably good results can be obtained even 
with a small SZP basis, particularly for lattice constants. 
This shows that structure optimizations with the LCAO 
code are likely to yield very accurate geometries. This 
is probably due to the fact that calculations of equilib- 
rium structures only involve energy differences between 
very similar structures, i.e. not with respect to isolated 
atoms, thus leading to larger error cancellations. 

With DZP the primary source of error in cohesive en- 
ergy comes from the free-atom calculation, where the 
confinement of each orbital raises the energy levels by 
around 0.1 eV. Thus, atomic energies are systematically 
overestimated, leading to stronger binding. This error 
can be controlled by using larger basis set cutoffs, i.e. 
choosing smaller orbital energy shifts during basis gener- 
ation. 



C. Structure optimizations 

LCAO calculations tend to reproduce geometries of 
grid-based calculations very accurately. In structure op- 
timizations, the LCAO code can therefore be used to pro- 
vide a high-quality initial guess for a grid-based structure 
optimization. 

While it is trivial to reuse a geometry obtained in one 
code for a more accurate optimization in another, our ap- 
proach is practical because the two representations share 
the exact same framework. Thus the procedure is seam- 
less as well as numerically consistent, in the sense that 
most of the operations are carried out using the same ap- 
proximations, finite-difference stencils and so on. With 
quasi-Newton methods, the estimate of the Hessian ma- 
trix generated during the LCAO optimization can be 
reused as well. For most non-trivial systems, an LCAO 
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FIG. 3: (Color online) PBE atomization energies from the G2-1 dataset, relative to the grid values. The corresponding Mean 
Absolute Errors with respect to the grid values are: 1.71 eV (20.4%) for Double Zeta (DZ); 0.36eV (4.45%) for Double Zeta 
Polarized (DZP); 0.25 eV (3.02%) for Triple Zeta Polarized (TZP)and 0.20 eV (2.44%) for Triple Zeta Double Polarized (TZDP). 



calculation is between 25 and 30 times faster than a grid 
calculation, making the cost of the LCAO optimization 
negligible. 

Fig. [5] shows a performance comparison when reusing 
the positions and Hessian from a LCAO-based structure 
optimization for a grid-based one, using the default ba- 
sis set. The system is a 38-atom truncated octahedral 
gold cluster with CO adsorbed, with the initial and final 
geometries shown in the inset. 

A purely grid-based optimization takes 223 CPU hours 
while a purely LCAO-based one, requiring roughly the 
same number of steps, takes 8.4 CPU hours. A further 
grid-based optimization takes 45 CPU hours, for a total 
speedup factor of 4. The value of an initial LCAO opti- 
mization is of course higher if the initial guess is worse. 
For systems where a large fraction of the time is spent 
close to the converged geometry, the speedup may not be 
as significant. 

The energy reference corresponds to the separate clus- 
ter and molecule at optimized geometries - the to- 
tal energy difference between an LCAO and a grid 



calculation is otherwise around 30 eV. It is therefore 
important to choose an optimization algorithm which 
will handle such a shift well. The present plots use 
the L-BFGS algorithm^^ (limited memory Broyden- 
Fletcher-Goldfarb-Shanno) from the Atomic Simulation 
Environment 



VI. CONCLUSIONS 

We have described the implementation of a localized 
basis in the grid based PAW code GPAW, and tested 
the method on a variety of molecules and solids. The 
results for atomization energies, cohesive energies, lat- 
tice parameters and bulk moduli were shown to converge 
towards the grid results as the size of the LCAO ba- 
sis was increased. Structural properties were found to 
be particularly accurate with the LCAO basis. It has 
been demonstrated how the LCAO basis can be used to 
produce accurate initial guesses (both for the electron 
wave functions, atomic structure, and Hessian matrix) for 



9 



o2 



Q 
I 

O 
< 
w 



0.3 



a (A) 





sz 




SZP 




DZ 




DZP 






SZ 


SZP 


DZ 


DZP 


GRID 


LiF 


4.08 


4.08 


4.02 


4.10 


4.06 


C 


3.61 


3.58 


3.59 


3.58 


3.57 


Na 


4.18 


4.19 


4.26 


4.24 


4.19 


MgO 


4.26 


4.28 


4.27 


4.27 


4.26 


Al 


4.24 


4.07 


4.08 


4.07 


4.04 


NaCl 


5.52 


5.62 


5.61 


5.67 


5.69 


Li 


3.68 


3.47 


3.70 


3.43 


3.43 


SiC 


4.50 


4.42 


4.46 


4.41 


4.39 


Si 


5.60 


5.52 


5.58 


5.49 


5.48 


A1P 


5.62 


5.55 


5.56 


5.53 


5.51 


Fc 


2.80 


2.77 


2.78 


2.83 


2.84 


Cu 


3.80 


3.59 


3.58 


3.64 


3.65 


PI 


4.02 


3.99 


3.95 


3.98 


3.98 



MAE 
MAE % 



0.097 
2.33 



0.034 
0.84 



0.068 
1.70 



0.019 
0.45 



0.5 



O.Ot 



I 

o 
< 

o 

^ -0.5H 



-1.0 



■ .1 .1 



I ■ I 



<-> 5 o < □ 
2 cn ^ in 



y 
55 



In i 
< 



u 



E a (cV) 





SZ 


SZP 


DZ 


DZP 


GRID 


LiF 


3.49 


4.48 


4.99 


4.52 


4.24 


C 


7.29 


7.51 


7.70 


7.89 


7.72 


Na 


0.97 


1.02 


1.07 


1.07 


1.09 


MgO 


2.81 


4.01 


4.94 


4.97 


4.95 


Al 


3.07 


3.51 


3.38 


3.54 


3.43 


NaCl 


2.94 


3.14 


3.24 


3.26 


3.10 


Li 


1.13 


1.58 


1.31 


1.63 


1.62 


SiC 


5.80 


6.31 


6.08 


6.48 


6.38 


Si 


4.14 


4.52 


4.34 


4.71 


4.55 


A1P 


3.77 


4.09 


3.92 


4.21 


4.08 


Fe 


1.34 


3.83 


4.77 


5.07 


4.85 


Cu 


2.38 


3.97 


3.75 


4.14 


3.51 


Pt 


4.54 


5.33 


5.57 


5.69 


5.35 


MAE 


0.86 


0.25 


0.19 


0.18 




MAE % 


20.70 


5.86 


5.51 


4.40 





in 



o < 



(J 



y 55 9: 
in < 



u 



B (GPa 





SZ 


SZP 


DZ 


DZP 


GRID 


LiF 


87 


84 


91 


70 


80 


C 


394 


408 


411 


422 


433 


Na 


8.9 


9.1 


8.3 


7.9 


7.9 


MgO 


156 


184 


209 


173 


154 


Al 


53 


74 


73 


79 


77 


NaCl 


35 


32 


34 


20 


24 


Li 


10.8 


15.2 


10.7 


16.3 


14.2 


SiC 


178 


196 


221 


202 


211 


Si 


70 


81 


77 


86 


88 


A1P 


69 


77 


76 


81 


82 


Fc 


248 


379 


297 


231 


198 


Cu 


88 


181 


166 


143 


141 


Pt 


224 


266 


309 


263 


266 


MAE 


22.9 


24.8 


23.2 


7.4 




MAE % 


20.4 


18.2 


18.8 


6.3 





FIG. 4: (Color online) Deviations in cohesive energy (top), lattice parameter (middle) and relative bulk modulus (bottom) 
from the converged grid results. The largest bars have been truncated and are shown with dotted edges - see the corresponding 
tables for the precise values. 
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Appendix A: Force formula 
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FIG. 5: (Color online) The energy as a function of iteration 
count (top) as well as CPU time (bottom) in structure opti- 
mizations. Shows a grid-based and an LCAO based structure 
optimization plus the continuation of the LCAO optimization 
after switching to the grid representation. 



subsequent grid-based calculations to increase efficiency 
of high-accuracy grid calculations. 

The combination of the grid-based and LCAO methods 
in one code provides a flexible, simple and smooth way 
to switch between the two representations. Furthermore 
the PAW formalism itself presents significant advantages: 
it is an all-electron method, which eliminates pseudopo- 
tential errors, and it allows the use of coarser grids than 
norm-conserving pseudopotentials, which increases effi- 
ciency. 

Finally, the LCAO method enables GPAW to perform 
calculations involving Green's function, which intrinsi- 
cally need a basis set with finite support. Current de- 
velopments along these lines include electron transport 
calculations, electron-phonon coupling and STM simula- 
tions. 



The force on atom a is found by taking the derivative 
of the total energy with respect to the atomic position 
R a . We shall use the chain rule on Eq. (|17j) . taking p^i/, 
-D?-, n(r), p(r), T^ v and v(r) to be separate variables for 
the purposes of partial derivatives: 



dE 



E 



+ E 



dE dp Vfl 
dp Vfi dR a 

SE 8n(r) 
6h(r) <9R Q 
dE 8T, 



v dE dD% 

DIJ J L 



dr 



SE dp(r) 



dT^ dR a 



Sp{r) dR a 
5E dv(r) 



dr 



Sv{r) dK a 



dr, (Al) 



where v(r) = ^ a v a (\r — R a |). The remaining quantities 
in the energy expression pertain to isolated atoms, and 
thus do not depend on atomic positions. The first term 
of Eq. (|Al~j) is 

^ d Pvil dR a ^ u ^ CmIn w 



dc 



y y QY^a S [iv c vn e nfm (A2) 



where we have used Eqs. (|29|) and (|36| in the first step, 
and Eq. (HIT) in the second. When the atoms are dis- 
placed (innnitesimally) , the coefficients must change to 
accommodate the orthogonality criterion. This can be 
incorporated by requiring the derivatives of each side of 
Eq. (|39[) to be equal, implying the relationship 



E V 



dS, 



dn a 



8r* 



\1V ILV 

Inserting this into Eq. (|A2|) yields 
\ - dE dp vll _ 



(A3) 



e Pvil dn a 



dn a 



where we have introduced the matrix 



(A4) 



(A5) 
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The equivalence of these forms follows from Eq. (|41[) . 
The overlap matrix elements depend on R a through 
the two-center integrals O^ and P^. The derivative of a 
two-center integral can be nonzero only if exactly one of 
the two involved atoms is a, and for nonzero derivatives, 
the sign changes if the indices are swapped. Taking these 
issues into account, Eq. (|A4I) is split into those three 
terms in Eq. (|44[) which contain E v)1 . 
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In the second term in Eq. (|A1|) . we take the 
dependent derivative for fixed p„ M , which by Eq. 
evaluates to 



D 



j 1 



dP b 



(A6) 



Again most of the two-center integral derivatives are 
A complete reduction yields the two terms in Eq. 
which depend on the vectors. 

Using Eq. ([50]). the third term of Eq. ([AT]) is 

SE dh(r) 



zero. 



Sn(r) dH a 



dr = 



v(t) „_,„ dr 



//// 

+ / u(r) 



9R 

<%"(|r 



<9R a 



■ dr. (A7) 



The sum over p can be restricted to p G a. 

Consider the fourth term of Eq. (IA1|) . Aside from 
n(r) and D^-, which are considered fixed as per the chain 
rule, the pseudo charge density p(r) depends only on the 



locations of the compensation charge expansion functions 
<7£(r) which move rigidly with the atom, so 



SE dp(r) 
Sp(r) dIL a 



dr 



5p(r) ^ SZ(r) d~g b L (r) 



X - uz 

SZ(r)j^W L (r) dR a 



dr 



a d~9l(r) 
dR a 



dr. 



(A8) 



The kinetic term from Eq. (|Aip is 
dE dT^ 



E 

/ 1 ' ' 



dTnu dR a 



(A9) 



and can also be restricted to p G a. Finally, the contri- 
bution from the local potential v a (r) is simply 



6E dv{r) 
5v(r) dK a 



dr = 



dv a (r - R a ) 
Hr) [ dR a j dr. (A10) 



By now we have considered all position-dependent vari- 
ables in the energy expression, and have obtained expres- 
sions for all terms present in Eq. 
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